Transition to Localization of Biased Walkers in a Randomly Absorbing Environment 
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We study biased random walkers on lattices with randomly dispersed static traps in one, two 
and three dimensions. As the external bias is increased from zero the system undergoes a phase 
transition, most clearly manifested in the asymptotic drift velocity of survivors which is zero below a 
critical bias. This transition is second-order in one dimension but of first order in higher dimensions. 
The model can be mapped to a stretched polymer with attractive interaction between monomers, 
and this phase transition would then describe sudden unfolding of the polymer when the stretching 
force exceeds a critical value. We also present precise simulations of the zero bias case where we 
show unambiguously that the transition between the Rosenstock and Donsker-Varadhan regimes is 
first order in dimension > 2. 
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I. INTRODUCTION 

The behavior of a particle diffusing amidst randomly 
dispersed static traps is already interesting by itself, but 
also because of its connections to several other problems: 
It is closely related to the Anderson model of delocaliza- 
tion and to some models of population dynamics [y . It 
can be mapped exactly onto flux lines in a medium with 
parallel columnar defects |]1| and to the collapse of a self- 
attracting polymer without excluded volume interaction. 

The problem is usually formulated as a random walk 
on a d— dimensional lattice; a random fraction c of all 
lattice sites is occupied by traps which kill a walker if it 
happens to land on them. The main quantity of interest 
is the survival probability after t steps, P{c,t), averaged 
over all walks and all trap configurations. For small t 
and small trap concentrations, P{c, t) decays exponen- 
tially as predicted by Rosenstock (RS) [||, but at later 
times it crosses over to the exact Donsker-Varadhan (DV) 
form |-| 



In P(c,i) 



^^^d+2id+2 



(1) 



where A = — ln(l — c) and Ad are exactly known nu- 
merical constants. Eq.(yJ) can be understood in terms of 
contributions from large trap-free regions. Such regions 
are very rare, but survival is much enhanced in them. 
Thus, P{c, t) is dominated by regions whose size grows 
with t as R^ f^ ji). 

The precise form of this crossover has been a matter of 
some controversy p-H]. In particular, simulations have 
been notoriously difficult for d > 2. While deviations 
from the RS prediction are easily seen |1, it has proven 
difficult to access the validity of Eq.p) in simulations 
p-||. For d = 3, the situation has largely been clarified 
by simulations of Anlauf llG] which unfortunately were 
never published, and by analytical work by Nieuwen- 
huizen 0. Very recently, similar scaling was found in 
0,|2| to hold also for d = 2. 

In this paper we first present further simulations of this 
model where we correct some previously made claims and 



verify the prediction of M that the RS-DV crossover is 
akin to a first order phase transition at infinite tempera- 
ture. We then show results of the model with a superim- 
posed bias. In the latter, the jumping rates in the (fixed) 
bias direction are enhanced by a factor & > 1 and those in 
the opposite direction are suppressed by the same factor. 

This external bias leads to a phase transition in the 
asymptotic dynamics. For small bias the long time dy- 
namics is similar to the zero-bias case: the walkers stay 
in large trap-free regions and do not drift. Instead, their 
typical positions increase slowly as x ~ i? ~ ii/(rf+2)_ g^i; 
above a critical bias b^ the surviving walks drift with a 
finite velocity ||l^. These behaviors are separated by a 
sharp transition which is of second-order in d = 1 but of 
ffist order in d > 2. At 6 \ 6c the drift velocity tends to 
zero in d = 1 but reaches a non-zero value i>c in d > 2. 
Asymptotically the survival probability decays exponen- 
tially in both phases, but the decay rate is singular at 
the transition point: its slope is discontinuous in d > 2 
while in d = 1 the second derivative is discontinuous. 

The averaging over disorder can easily be done exactly 
and leads to 



with 



Here, 



P(c,i) = e-T'(2d)-*Z, 



= Y^ii-cYb^ 

walks 



7 = ln( 



b + b- 



2d 



1) 



(2) 



(3) 



(4) 



and s is the number of distinct sites visited. Finally, 
Zt can be understood as the partition sum of a self- 
attracting polymer (without excluded volume interaction; 
we assume ksT = 1) of length t, stretched by applying 
opposite forces ±ln6 to its two ends. The self attrac- 
tion comes about by the fact that each monomer which 
is placed on a new site is punished by a Boltzmann factor 
1 — c, while no such factor is applied when a monomer 



is placed at a site which is aheady occupied by another 
monomer. 

Although this is rather artificial as a polymer model, 
the mapping is very useful for suggesting efhcient simu- 
lation methods and for interpreting the results. 

Our simulations are performed using the Pruned- 



Enriched-Rosenbluth Method (PERM) [Ij which is a 
growth scheme to build weighted walks: Walks with too 
small weights are pruned, promising configurations are 
further encouraged by cloning them if their weight ex- 
ceeds a (length-dependent) given value. PERM has been 
used to simulate efficiently a large number of polymer 
problems including 0-collapse, DNA denaturation, pro- 
tein folding, and critical unmixing 1 15|. In addition to the 
"population control" provided by cloning and pruning, 
walks are also guided by assigning different a priori prob- 
abilities to different jumps. In particular, we chose these 
such that previously unvisited sites are unfavored against 
previously visited ones, Punvisited = (1 - c)pmsited, and 
jumps in the bias direction are more likely than others, 
p+x '■ P-x '■ P± = Vb : ■\/l/& : 1. These a priori biases 
are of course compensated by appropriate weight factors 
oc l/p. 



II. UNBIASED WALKERS 

For b = 1 walks are asymptotically subdiffusive or lo- 
calized; the DV theory predicts the mean squared dis- 
placement to be (i?^) ~ i2/(d+2)^ However, for small 
trap concentrations walks are brownian at small times 
((i?^) ^ t) reverting to DV scaling asymptotically from 
above [ |lO| . This curious nonmonotonici ty (als o observed 
in other models of polymer collapse ||l^,|lj,0 results be- 
cause the most successful walks at shorter times do not 
survive for longer periods where DV scaling holds (Fig. 
la). The same is seen, albeit much less dramatically, in 
d = 2 (Fig. lb). It is not seen in d = 1 where walks 
remain localized for all times and (R'^) is monotonic. 
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FIG. 1. (a) (i?^) vs t in d = 3 for unbiased walks, showing 
the sudden localization indicated by the drop in {-R^). The 
time at which this localization occurs decreases with increas- 
ing trap concentration. The transition proceeds more gradu- 
ally for larger c, and for c >= 0.7 no maximum is. (b) Similar 
data for d — 2. The collapse is there much less pronounced. 



The sharp collapse in d = 3 suggests that the RS-DV 
crossover is indeed a first order phase transition in the 
limits i — > CX3 and zero trap concentration, as predicted 
in H. In the polymer model, it would be a first order 
transition at infinite temperature. A more direct sup- 
port of this is provided by histograms of the distribution 
of visited sites. For a first order transition we expect two 
peaks for large i, the positions of which are separated by 
a distance c>c t. This is indeed seen in d = 3 (Fig. 2), but 
not for d = 2. The latter is not surprising in view of the 
much softer collapse in d = 2. It could just mean that 
the two-peak structure develops only for t values beyond 
those obtainable by our simulations (i « 10^). 

Combining the RS and DV limits into a single scaling 
ansatz, Anlauf [Q suggested for d = 3 that the average 
number of visited sites scales as 



t-\s) = g,iXt'/') 



(5) 



A similar ansatz was made for d = 2 in ||ll| , but as shown 
in |12| the proper ansatz for d = 2 must involve logarith- 
mic corrections. 



t-i(s)lnt = 52(At^''Vlni) • 



(6) 



As seen from Fig. 3a, Eq.(0) does not hold with the sim- 



ple scaling function proposed in |10|. Anlauf 's data had 
suggested a piecewise linear 33(0;) on a log-log plot, with 
53 (x) — const for x < x* and gsix) = const/ x^'^ for 
X > X* . We see instead from Fig. 3a that gz{x) develops 
a kink with infinite slope at a; = x* . This is not too 
surprising. While a behavior as claimed in |l^ would 
be natural for a second order transition, a more singu- 
lar behavior is expected for a first order transition. For 
d = 2 the scaling function is much smoother (see Fig. 3b), 
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FIG. 2. Histogram of visited sites for unbiased walks in 
d = 3. The two peaks correspond to the collapsed (DV, left) 
and extended (RS, right) phases. Trap concentrations (resp. 
monomer interaction strenghts) have been choosen such that 
both peaks have equal heights; c = .01993 (t = 8192), .01283 
(16384), .00823 (32768), .00527 (65536). 



although this might be a transient effect only. We finally 
point out that we require the precise asymptotics for 2- 
d random walks as derived in [Q , with In t replaced by 
ln(8i), to obtain a good data collapse. 
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FIG. 3. (a) {s)/t vs At^/^ in d = 3 for unbiased walks. For 
large t the data collapse onto a scaling function which is flat 
for small arguments, oc (Ai^")~^'^ for large arguments, and 
has a cusp with infinite right slope in between. It is not clear 
from the data whether the scaling function is continuous or 
discontinuous at the cusp, (b) Similar for d = 2. Notice that 
we replaced Ini in Eq.(6) by ln(8t) in order to improve the 
data collapse. The straight line is the asymptotic DV predic- 
tion. 
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FIG. 4. (a) Scaling plot for P{c, t) for unbiased walks in 
d — 3. The straight lines are the RS and DV predictions. 
While the RS prediction is realized numerically (for small con- 
centrations c) up to the transition point, the true DV scaling 
is reached only very slowly, (b) Similar for d — 2. 



Using Eq.(H), Anlauf had proposed a similar scaling 



t-i/3lnP(c,t) = /3(Ai2/3). 



(7) 



Using his simple gs (x) he derived, by thermodynamic in- 
tegration, an /a (a;) which gave a very slow convergence to 
the asymptotic DV scaling after having deviated sharply 
and suddenly from RS behavior. Although we disagree 
with his scaling function, our results for P{c, t) (Fig. 4a) 
fully confirm with his predictions within the numerical 
errors. Notice that this disagrees slightly with the re- 
sults found in [Q . Again the analogous ansatz for d = 2 
should involve logarithms |12| , and for a faster conver- 
gence we must replace Int by ln(8t). 



lnP(c,i)/ln(8i) ^ f2{Xt^/^/H^t)) 



(8) 



Agreement with simulations is very good, and again we 
see very slow convergence to the DV asymptotics, in con- 
trast to ill. 



III. NONZERO BIAS 

For any nonzero bias the long time probability decay 
is exponential [^,^. This is most easily seen as fol- 
lows: Let (/)(a;,x_L, i; 6;C) be the probability density, for 
any fixed trap configuration C and bias b, that a walker 
who started at i = in the origin x = is at time t at 
the position x = (a;,x^). One verifies easily that this is 
related to the analogous distribution of the bias-free case 
by (see [El| for the continuum case) 



[X,li± 



i;6;C) = e"'^*6^0(x,Xi,i;l;C) 



(9) 



This 'gauge transformation' means in particular that the 
return probability 



Pretur„W = (</>(0,O,t;6))c 

decays exactly as 



Pr. 



.(i) 



ttp{b=l) 



it), 



(10) 



(11) 



i.e. as an exponential multiplied by the DV stretched ex- 
ponential since P4^~^„(i) also satisfies the DV estimate 
up to algebraic factors. 

For the survival probability we obtain 

In P{c, t) = In P(^=i) (c, t)--ft + ln(cosh(a; In 5)) walks, C • 

(12) 

Since the last term is positive definite, we have the exact 
inequality In P(c,i) > In p(^=^)(c,t) — ^t. Moreover, if 
the walkers stay localized in spite of the bias, x increases 
sublinearly with t and thus the asymptotic decay rate is 



a = — lim - In P{c, <) = 7 

t — 'OO t 



(localized) (13) 



Notice that this is independent of the trap concentra- 
tion. A smaller decay rate can be obtained only when 
the walkers drift, but then the chance to hit a trap is 
finite at any time step, and P{c,t) decays again expo- 
nentially. Indeed, in the drifting regime any correlations 
should be short, and thus both the decay rate and the 
drift velocity are attained already after a short time. 



predicted theoretically. Secondly, we see that all curves 
for b < 1.11 become horizontal for t —^ oo, in agreement 
with the prediction be — 1/(1 — c) = 1.111 . . .. The 
decay rates obtained from Fig. 5 for b > be are shown in 
Fig. 6a, together with the theoretical predictions. We see 
that the transition is indeed smooth, as obtained theo- 
retically. Finally, the drift velocities for c = 0.1 are shown 



A. d= 1 

In one dimension it is possible to derive analytic esti- 
mates for the survival probability and the drift velocity of 
delocalized walks. An exact lower bound on the survival 
probability in terms of probabilities of large trap-free re- 
gions is obtained in a way similar to B , by estimating the 
contributions of trap-free regions of length I downstream 
of the starting point of the walker. Take some positive 
integer I and consider configurations with no trap in the 
interval [0,1]. We then have 



p{c,t)>{i~cy Pi{t) 



(14) 



where the first factor on the rhs. is the probability to 
find such a gap in the trap configuration, and Pi{t) is the 
probability to survive in it. The latter can be estimated 
using Eq.(ffl, the method of images to incorporate the 
boundary conditions at a; = — 1 and x = I + 1, and Stir- 
ling's formula. We obtain, up to a constant dependent 
only on b, 



\nPi{t) > llnb-'^t- 






^n(l-i) 

(15) 



The rhs. ofEq.(^ 
c)6)] when b > 



|) is maximal iov I = I* =t tanh[ln((l — 
be = 1/(1 — c), and for / — o{t) 
when b < be- Thus the decay rate is bounded by 
a < 7 — In [cosh (In (1 — c)b))] for b > be- Actually we 
claim that this bound is saturated, since it can also 
be understood as the result of a (non-rigorous) sad- 
dle point approximation. Summarizing, we have thus 
found that the decay rate is equal to 7 for b < be and 
is given by a = 7 — ln[cosh(ln(l — c)b))] for b > be- 
Similarly, the drift velocity of the survivors, defined as 

V = lim4^oo(2;)/i, vanishes for b < be and is given by 

V = tanh[ln((l — c)b)] for b > be- In the continuum limit 
(c -> 0, {1- c)b = finite) this agrees with [|2|,|l|,|3| , 
but for finite c it disagrees with p2| [M. Thus there 
is a phase transition at 6 = 6c- It is second order in the 
sense that v, a, and da/ db are continuous aXb — be, while 
dv/db and d^a/db^ are discontinuous. 

These predictions are fully confirmed by simulations, 
as seen from Figs. 5 and 6. In Fig. 5 we see In P(c, t) + -ft 
versus t for c = 0.1 and for several values of 6. First of 
all we see that all curves are above that for 6 = 1, as 
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FIG. 5. In P(c,i) + 7t for d — 1 and trap concentration 
c — 0.1. Statistical errors are smaller than the thickness of 
the curves. 
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FIG. 6. (a) Decay rate a as a function of bias 6, for c = 0.1 
{d = 1). For h > 1.2 the statistical errors are much smaller 
than the size of the symbols, for smaller b there is an increas- 
ing systematic uncertainty due to the non-trivial extrapola- 
tion for t ~> oo- (b) Analogous results for the drift velocity. 
Dashed curves are analytic predictions. 



in Fig. 6a. Again agreement with theory is perfect. We 
should add that similar results were also found for other 
values of c, and that (x) increased always monotonically 
with t, in contrast to what is found for d > 1 (see below). 



d — 2 and c = 1/8) we see clearly two peaks in P{x), 
for t > 900. The left peak is at x/t « 0. Indeed, it is 
slightly to the right of x/t — 0, but it moves to x/t — 
for i — > oo. It of course corresponds to the localized 
phase. The right peak is the drifting phase. The fact 



B. d> 2 



In higher dimensions we don't have similar analytic 
results. Moreover, direct simulations of the transition re- 
gion are very difhcult, even with our efficient algorithm. 
But Eq. ( p^ ) together with the expected fast convergence 
of the decay rate and drift velocity in the delocalized 
phase allow us, in spite of this difhculty, to compute w, a, 
and be with very high precision. 

In Fig. 7 we show the average displacements (x) in 
d — 2 against t, for c = 1/8 and for several values of 
b. Very similar data are obtained for d = 3 and for other 
values of c. Obviously be > 1.65, since for all b < 1.65 
the curves decrease after an initial rise (x) oc t. The 
reason for this non-monotonic behavior is the same as in 
Fig.l: Walkers who venture far out do better initially, 
but finally only those win who started in a large trap- 
free region and stayed in it. From Fig. 7 one might guess 
that 1.65 < be < 1.7, but this would be wrong: Also 
for b = 1.7 the curve should bend down ultimately, but 
this will happen very late (see below) and our algorithm 
might easily miss it (in the transition region it generates 
many drifting walks with low weight and only very few 
compact walks, but these with high weight). 

The true transition point can be estimated from the 
survival probabilities, shown in Fig. 8 for d = 2. Actu- 
ally, in view of eq.(p^), we plot In P(c, t) + jt. For small 
b we see curves which follow closely the curve for 6=1, 
but are slightly above it. These curves contain a linear 
part at small t whose slope decreases and whose length 
increases with b. The cross-over from this linear part to 
the part following the DV curve is rather sharp and stays 
so with increasing b. Since the DV curve is asymptoti- 
cally flat, this scenario should hold up to the value of 
b where the "initial" part becomes horizontal. Beyond 
that value we are in the drifting phase. Obviously the 
initial straight parts in Fig. 8 correspond precisely to the 
(transient) drifting phase seen already in Fig. 7. Indeed, 
the cross-overs between the two regimes are at the same 
positions in Figs. 7 and 8. 

This means that be is implicitly given by the con- 
dition ttdrifting = 7, whcre adrifting IS the decay rate 
in the drifting phase. While this phase is asymptot- 
ically stable for Odrifting < 7, it is only transient for 
cidrifting > 7- Evcn in the transient regime it can be 
measured precisely (together with the drift velocity), if 
this regime lasts long enough. Moreover, the drifting 
phase exists also for b < be and large t, but is only 
metastable in this regime. This is clearly seen from his- 
tograms of the x-distribution. In Fig. 9 (which refers to 
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FIG. 7. Average displacement in the bias direction for 
c = 1/8 in d = 2. 
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FIG. 8. In P(c,i) + 7t for d — 2 and trap concentration 
c — 1/8. Since these curves must be increasing with b, the 
fact that the curve for b = 1.65 crosses for t > 6000 below 
those for 6 = 1.555 to 1.635 is due to numerical inaccuracy. 



that these peaks become sharper and move apart with 
increasing t shows that the transition is first order. For 
d > 2 we expect the transition then to be of first order a 
fortiori. 

Since 7 is known analytically, we can locate be by sim- 
ulations in the delocalized phase alone. This procedure 
is indicated in Fig. 10 where the dashed curve is 7(6) and 
the various continuous curves arc finite-time decay rates 
obtained from short-time runs with different trap concen- 
trations. The error bars of the latter are much smaller 
than the size of the symbols. The intersection points give 
be and ac in terms of the trap concentration c. 



The scaling behavior for c ^- is easily obtained by 
observing that 7 « (6 — l)^/(2(i) for 5 « 1, while a ^ c 
in the drifting phase (see Fig. 11). The latter is true for 
small c since the Rosenstock approximation becomes ex- 
act in this limit. Therefore we predict that be — 1 ^ \fc. 
This is verified in Fig. 12 for d = 2 and d = 3. Finally, 
from Fig. 13 we see that also v^. ~ 'Jc, i.e. Vc oc h^. The 
latter is very natural, but we must remember that it is 
not true in d = 1 . 
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FIG. 9. Histograms of the rescaled end point distance 
P{x/t) versus x/t for d = 2 and c = 0.125. Bia- 
ses were adjusted so that both peaks have equal height: 
h = 1.510 (i = 900), 1.577 (i = 1800), 1.608 {t = 2700), 1.628 
(i — 3600). Normalization is arbitrary. The peak at x/t ~ 
is due to the localized phase, the other corresponds to the 
drifting phase. 
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FIG. 11. Log- log plot of the critical decay rate versus c. 
For small c, both curves scale oc c. 
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FIG. 12. Log-log plot of the logarithm of the critical bias 
versus c. For small c, both curves scale as ^/c. 
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FIG. 10. Decay rates a measured in the drifting (delo- 
calized) phase, d = 3. The dashed curve is the prediction 
a = 7(6) for the localized phase. To the left of this curve, the 
drifting phase is only transient. In spite of that, a is easily 
measured with high precision, thus the statistical errors are 
< 10"'*, much smaller than the symbols. 
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FIG. 13. Log-log plot of the critical drift velocity versus c. 



IV. ALTERNATIVE INTERPRETATIONS 
A. Population Dynamics 

Up to now, we dealt with independent walkers. De- 
tailed studies of this model were often criticized on the 
grounds that the tiny probabilites that obtain in the DV 
regime render the results physically irrelevant (in the 
model with bias, probabilities are even smaller). This 
criticism can be countered by adding autocatalytic par- 
ticle production (local production of off-springs) |^ . 

Thus we consider a model where particles can hop, get 
absorbed when hitting a trap, and can reproduce with 
some fixed rate a. For simplicity we do not include any 
interaction between particles, in particular we do not take 
into account any excluded volume effect or any compe- 
tition for food or other resources. The only effect of re- 
production is then that the average number of survivors 
at generation t is P{c,t){l + aY . 

In the bias-free case any non-zero a will thus lead to an 
explosion of the population, indicating that the model is 
basically sick. But in the case with drift a suitably cho- 
sen rate a will lead to either a much reduced decay of 
the particle number, or even to a stationary population 
size. In this case the neglect of excluded volume and sat- 
uration effects will be much more benign, and the effects 
studied in this paper might have a bigger chance to be 
observable in real experiments. 

Finally, including excluded volume effects would lead 
to directed percolation (contact process) in a random 
medium |E5|. 



the collapse of the unstretched polymer happens at infi- 
nite temperature in the present model. In the standard 
polymer models with next-nearest neighbor attraction or 
with three body repulsion the O-coUapse is second order, 
in contrast to the present model. But the dependence 
of (i?^) on chain length is also non-monotonic in d > 3 
p6| , p^jT7[. In d = 4 the collapse is indeed "quasi first 
order" |17|] , and special care is needed to realize that it 
actually is continuous. We believe that in the present 
model the transition is truely of first order. 



C. Magnetic Flux Lines 

A last interpretation is in terms of magnetic flux lines 
in a type II superconductor with parallel columnar de- 
fects ]l|. These defects are quenched and randomly dis- 
tributed. They tend to pin the flux lines, thus preserving 
superconductivity in the sample. In the unbiased case, 
the mean orientation of the flux lines is parallel to the 
defects, caused by a magnetic field in the same direction. 
The biased case corresponds then to a non-zero angle be- 
tween the field and the defects. For small magnetic field 
the flux lines will essentially follow the defects. But above 
a critical field strength (resp. tilting angle) the average 
angle between flux lines and defects is non-zero. Our 
main result is that this transition is flrst order in the 
sense that the average angle jumps by a flnite amount 
when passing through the critical field strength. This 
jump scales as the square root of the defect density. 



B. Stretched Collapsed Polymers 



V. CONCLUSION 



We have already pointed out that Zt (Eq.(p|)) can be 
understood as the partition sum of a polymer with self- 
attraction of strength ksT X, and stretched by a force 
ksT Inb. The decay rate is then essentially the free en- 
ergy per monomer in the thermodynamic limit. 

The fact that unfolding of a collapsed polymer is a first 
order transition is well known pq-p8[. The coexistence 
of stretched and coiled phases corresponds to a "tadpole" 
configuration py,M where part of the chain is collapsed 
and the rest sticks out of the coil, so that the total length 
X is equal to some prescribed value. Our finding that the 
critical bias is ex ^/c means then that the string tension 
in a tadpole configuration scales, in the limit of a very 
long chain and close to the collapse point, as the square 
root of the monomer-monomer attraction strength. This 
seems to be a new result, and it remains to see whether 
it holds also for more realistic polymer models. 

The main difference to standard polymer models is the 
fact that we have neglected here excluded volume effects 
(they have no counterpart in the trap problem). Thus 



Using extensive Monte Carlo simulations by means of 
the PERM algorithm, we have studied the problem of 
particles performing unbiased and biased random walks 
in media with randomly located traps. In the bias free 
case we find the first clear numeric indication for the 
first order nature of the Rosenstock - Donsker-Varadhan 
transition in 3 dimensions. In d = 2 the situation is 
less clear, but we showed that also there the transition 
is rather sharp. When the diffusion is biased, we find a 
transition from localized to drifting walkers at a critical 
bias which depends on the trap density. This transition 
is second order in d = 1 but first order in d > 2. Com- 
bining the simulations with exact analytic estimates we 
were able to draw the phase diagram with high accuracy. 
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